#install.packages("quantreg")

directory<-Sys.getenv("US_Ineq_Repl")
MMwd<-file.path(directory, "Scripts/MM")
Datawd<-file.path(directory, "Processed")
Figwd<-file.path(directory, "Results/Figures")
  
setwd(Datawd)
load("QR_mod.RData")

require(quantreg)
setwd(Figwd)

# union
ii<-which(rownames(model0$coefficients)=="union")
jj<-which(rownames(model1$coefficients)=="union")

png('QR_Coef_UnionB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Unionization", 
     ylim=c(-0.05,0.35),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
mtext("Quantile",side=1,line=1)

par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="",
     ylim=c(-0.05,0.35),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), 
     pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)

legend(0.05,0.1,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

# Public Sector
ii<-which(rownames(model0$coefficients)=="pubsect")
jj<-which(rownames(model1$coefficients)=="pubsect")

png('QR_Coef_PublicSectorB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Public Sector", 
     ylim=c(-0.10,0.10),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
mtext("Quantile",side=1,line=1)

par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="",
     ylim=c(-0.10,0.10),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), 
     pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)

legend(0.05,0.1,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()


# Manufacturing Sector
ii<-which(rownames(model0$coefficients)=="manuf")
jj<-which(rownames(model1$coefficients)=="manuf")

png('QR_Coef_ManufacturingSectorB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Manufacturing Sector", 
     ylim=c(-0.10,0.45),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
mtext("Quantile",side=1,line=1)

par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="",
     ylim=c(-0.10,0.45),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), 
     pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)

legend(0.05,0.1,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

#nonwhite
ii<-which(rownames(model0$coefficients)=="nonwhite")
jj<-which(rownames(model1$coefficients)=="nonwhite")

png('QR_Coef_NonWhiteB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Nonwhite", 
     ylim=c(-0.2,0),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(-0.2,0),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)

mtext("Quantile",side=1,line=1)
legend(0.2,-0.02,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()


#female
ii<-which(rownames(model0$coefficients)=="female")
jj<-which(rownames(model1$coefficients)=="female")

png('QR_Coef_FemaleB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Women", 
     ylim=c(-0.3,-0.05),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(-0.3,-0.05),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)
mtext("Quantile",side=1,line=1)
legend(0.3,-0.1,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

# Part time
ii<-which(rownames(model0$coefficients)=="partte")
jj<-which(rownames(model1$coefficients)=="partte")

png('QR_Coef_PartTimeB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Part Time", 
     ylim=c(-0.3,-0.05),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(-0.3,-0.05),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)
mtext("Quantile",side=1,line=1)
legend(0.3,-0.1,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

#married 
ii<-which(rownames(model0$coefficients)=="married")
jj<-which(rownames(model1$coefficients)=="married")

png('QR_Coef_MarriedB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Married", 
     ylim=c(0,0.12),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,0.12),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)
mtext("Quantile",side=1,line=1)
legend(0.3,0.04,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

# Urban
ii<-which(rownames(model0$coefficients)=="smsa")
jj<-which(rownames(model1$coefficients)=="smsa")

png('QR_Coef_UrbanB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Urban", 
     ylim=c(0,0.15),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,0.15),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)
mtext("Quantile",side=1,line=1)
legend(0.3,0.04,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

# experience
ii<-which(rownames(model0$coefficients)=="exper")
jj<-which(rownames(model1$coefficients)=="exper")

png('QR_Coef_ExperienceB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Experience", 
     ylim=c(0.0,0.045),col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,0.045),col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)  
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)
mtext("Quantile",side=1,line=1)
legend(0.1,0.04,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()


# High School
ii<-which(rownames(model0$coefficients)=="educ2")
jj<-which(rownames(model1$coefficients)=="educ2")

png('QR_Coef_HighSchoolB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="High School", 
     ylim=c(0,1.2) ,col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,1.2) ,col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)  
mtext("Quantile",side=1,line=1)
legend(0.1,0.5,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()


# Some College
ii<-which(rownames(model0$coefficients)=="educ3")
jj<-which(rownames(model1$coefficients)=="educ3")

png('QR_Coef_SomeCollegeB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Some College", 
     ylim=c(0,1.2) ,col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,1.2) ,col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)  
mtext("Quantile",side=1,line=1)
legend(0.1,0.7,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

# Associate Degree
ii<-which(rownames(model0$coefficients)=="educ4")
jj<-which(rownames(model1$coefficients)=="educ4")

png('QR_Coef_AssociateDegreeB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="Associate Degree", 
     ylim=c(0,1.2) ,col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,1.2) ,col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)  
mtext("Quantile",side=1,line=1)
legend(0.1,0.7,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()


# College or more
ii<-which(rownames(model0$coefficients)=="educ5")
jj<-which(rownames(model1$coefficients)=="educ5")

png('QR_Coef_CollegeB.png',width = 500, height = 500)
plot(smodel0, parm = ii, ols=F, xlab = "",ylab = "",main="College Degree", 
     ylim=c(0,1.2) ,col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
par(new=TRUE)
plot(smodel1, parm = jj, ols=F, xlab = "",ylab = "",main="", 
     ylim=c(0,1.2) ,col=c(rgb(100,150,190,70,maxColorValue=255),rgb(100,150,190,70,maxColorValue=255)), pch = 19, cex=2,
     type="l",lwd = 2)
lines(model1$tau,model1$coef[jj,],lty="dashed",col="blue2",lwd=2)  
mtext("Quantile",side=1,line=1)
legend(0.2,0.4,
       legend=c(expression(hat(beta)[1986] (tau)) ,"95% C.I.",expression(hat(beta)[2015] (tau)) ,"95% C.I."), 
       bty = "n",
       col=c("chocolate1",rgb(230,130,190,80,maxColorValue=255),"blue2",rgb(100,150,190,70,maxColorValue=255)), 
       lwd =c(2,8,2,8), 
       lty=c("solid","solid","dashed","solid"), 
       x.intersp=c(1,1,1,1),
       ncol = 2, 
       cex = 1.1
)
dev.off()

rm(list = ls())
